Adjusting the Intercept

One of the ways of adjusting your logistic regression model to account for the sampling to balance the events with non-events is by adjusting the intercept from the logistic regression model. The intercept is what sets the “baseline” probability in your model which is now too high. You can adjust your predicted probabilities from the regression model to account for this incorrect intercept. They can be adjusted with the following equation:

\[ \hat{p}_i = \frac{\hat{p}_i^* \rho_0 \pi_1}{(1-\hat{p}_i^*)\rho_1\pi_0 + \hat{p}_i^* \rho_0 \pi_1}\]

where \(\hat{p}_i^*\) is the unadjusted predicted value, \(\pi_1\) and \(\pi_0\) are the population proportions of your target variable categories (1 and 0), and \(\rho_1\) and \(\rho_0\) are the sample proportions of your target variable categories (1 and 0).

Luckily SAS will do this for us! All we have to do is give SAS the prior proportions from our population. In PROC FREQ we are again looking at the original distribution of the churn variable, but this time recording the proportions of 1’s and 0’s. The out = option saves the frequency analysis to a data set here called priors. In this data set we drop the percent variable with the drop = option and rename the count variable to prior with the ```rename = ```` option.

From there we build our logistic regression in PROC LOGISTIC as we have previously learned in earlier sections. The main difference comes in the SCORE statement. The data = option specifies the data set we want to score (here the validation data set) and the out = option specifies the name of our newly scored data set. The prior = option specifies the name of the data set that contains the prior proportions for the target variable (the priors data set we created from PROC FREQ).

proc freq data=logistic.tele_churn noprint;
   table churn / out=priors(drop=percent rename=(count=_prior_));
run;
        
proc logistic data=churn_t;
   class international_plan(ref='no') voice_mail_plan(ref='no') / param=ref;
   model churn(event='TRUE') = international_plan
                               voice_mail_plan 
                               total_day_charge 
                               customer_service_calls / clodds=pl;
   score data=churn_v prior=priors out=churn_scored1;
run;
quit;

This will produce predicted probabilities for our validation data set than have been adjusted to account for the sampling we did to balance out the data set.

Which approach is better? The more common approach is the weighted observation approach. Empirical simulation studies have proven that for large sample sizes (n > 1000), the weighted approach is better. In small sample sizes (n < 1000), the adjusted intercept approach is only better when you correctly specify the model. In other words, if you plan on doing variable selection because you are unsure if you have the correct variables in your model, then your model may be misspecified in its current form until after your perform variable selection. That being the case, it is probably safer to use the weighted approach. This is why most people use the weighted approach.

Weighted Observations

The other technique for adjusting the model for over or under sampling is by building the model with weighted observations so the adjustment happens while building instead of after the fact. In weighting the observations we use weighted MLE instead of plain MLE since each observation doesn’t have the same weight in the estimation of the parameters for our model. The only question now is what are the weights. We derive the weights from the same formula we had for the intercept adjustment. The weight for the rare event is 1, while the weight for the non-event is \(\rho_1\pi_0/\rho_0\pi_1\). For our data set this makes the weights 1 and 18.49 respectively.

We use the DATA STEP to create a weights variable in our training data set with these weights. Next, we build the logistic regression model for our training data as we previously learned in the last section, but with the new WEIGHT statement. In this WEIGHT statement we specify the weights variable for PROC LOGISTIC to use these weights for the weighted MLE.

data churn_t;
   set churn_t;
   weights = 1;
   if churn = 'FALSE' then weights = 18.49;
run;

proc logistic data=churn_t;
   class international_plan(ref='no') voice_mail_plan(ref='no') / param=ref;
   model churn(event='TRUE') = international_plan
                               voice_mail_plan 
                               total_day_charge 
                               customer_service_calls / clodds=pl;
   weight weights;
   score data=churn_v out=churn_scored2;
run;
quit;

Now we can analyze our results above knowing that the predictions and model have already been adjusted for our sampling because of the weighted approach.

Which approach is better? The more common approach is the weighted observation approach. Empirical simulation studies have proven that for large sample sizes (n > 1000), the weighted approach is better. In small sample sizes (n < 1000), the adjusted intercept approach is only better when you correctly specify the model. In other words, if you plan on doing variable selection because you are unsure if you have the correct variables in your model, then your model may be misspecified in its current form until after your perform variable selection. That being the case, it is probably safer to use the weighted approach. This is why most people use the weighted approach.

Adjusting the Intercept

One of the ways of adjusting your logistic regression model to account for the sampling to balance the events with non-events is by adjusting the intercept from the logistic regression model. The intercept is what sets the “baseline” probability in your model which is now too high. You can adjust your predicted probabilities from the regression model to account for this incorrect intercept. They can be adjusted with the following equation:

\[ \hat{p}_i = \frac{\hat{p}_i^* \rho_0 \pi_1}{(1-\hat{p}_i^*)\rho_1\pi_0 + \hat{p}_i^* \rho_0 \pi_1}\]

where \(\hat{p}_i^*\) is the unadjusted predicted value, \(\pi_1\) and \(\pi_0\) are the population proportions of your target variable categories (1 and 0), and \(\rho_1\) and \(\rho_0\) are the sample proportions of your target variable categories (1 and 0).

Luckily, this is rather easy to do in R as we can directly adjust the predicted probability vector with this equation. First, we build a logistic regression model as seen in the previous section using the glm function on our undersampled training data. We can use the predict function to obtain our predicted probabilities from this model. These are biased probabilities however because of our sampling. By simply applying the equation above on the predicted probability vector we can get unbiased predicted probabilities. We then store these as a new column in our validation data set.

logit.model <- glm(churn ~ factor(international.plan) + 
                           factor(voice.mail.plan) + 
                           total.day.charge + 
                           customer.service.calls, 
                   data = train_u, family = binomial(link = "logit"))
summary(logit.model)

Call:
glm(formula = churn ~ factor(international.plan) + factor(voice.mail.plan) + 
    total.day.charge + customer.service.calls, family = binomial(link = "logit"), 
    data = train_u)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0049  -0.7499  -0.0596   0.8002   2.1917  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -4.71202    0.77589  -6.073 1.25e-09 ***
factor(international.plan)yes  2.91300    0.58964   4.940 7.80e-07 ***
factor(voice.mail.plan)yes    -0.30174    0.43242  -0.698    0.485    
total.day.charge               0.09769    0.01829   5.341 9.26e-08 ***
customer.service.calls         0.63437    0.12268   5.171 2.33e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 296.67  on 213  degrees of freedom
Residual deviance: 217.88  on 209  degrees of freedom
AIC: 227.88

Number of Fisher Scoring iterations: 4

valid_p_bias <- predict(logit.model, newdata = valid, type = "response")
valid_p <- (valid_p_bias*0.5*(154/3004))/((1-valid_p_bias)*0.5*(2850/3004)+valid_p_bias*0.5*(154/3004))

valid <- data.frame(valid, 'Pred' = valid_p)

head(valid_p)

          2           5           9          11          16          17 
0.009815246 0.488633041 0.019242082 0.010848222 0.071763867 0.015565721 

We then store these as a new column in our validation data set.

Weighted Observations

The other technique for adjusting the model for over or under sampling is by building the model with weighted observations so the adjustment happens while building instead of after the fact. In weighting the observations we use weighted MLE instead of plain MLE since each observation doesn’t have the same weight in the estimation of the parameters for our model. The only question now is what are the weights. We derive the weights from the same formula we had for the intercept adjustment. The weight for the rare event is 1, while the weight for the non-event is \(\rho_1\pi_0/\rho_0\pi_1\). For our data set this makes the weights 1 and 18.49 respectively.

This is easily done in R by creating a new weight variable in our training data with the ifelse function. If our churn variable takes a value of TRUE then that observation has a weight of 1, while the other observations have a weight of 18.49. From there we build the logistic regression model as we learned in the previous section using the glm function. The only new option is the weights = option where you specify the vector of weights for the training data set so weighted MLE can be performed.

train_u$weight <- ifelse(train_u$churn == 'TRUE', 1, 18.49)

logit.model.w <- glm(churn ~ factor(international.plan) + 
                             factor(voice.mail.plan) + 
                             total.day.charge + 
                             customer.service.calls, 
                     data = train_u, family = binomial(link = "logit"), 
                     weights = weight)
summary(logit.model.w)

Call:
glm(formula = churn ~ factor(international.plan) + factor(voice.mail.plan) + 
    total.day.charge + customer.service.calls, family = binomial(link = "logit"), 
    data = train_u, weights = weight)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.6377  -0.9171   0.3496   2.0589   3.3592  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -7.29928    0.52017 -14.033  < 2e-16 ***
factor(international.plan)yes  2.57544    0.28830   8.933  < 2e-16 ***
factor(voice.mail.plan)yes    -1.08471    0.33129  -3.274  0.00106 ** 
total.day.charge               0.10566    0.01363   7.754 8.90e-15 ***
customer.service.calls         0.42969    0.07811   5.501 3.78e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 843.97  on 213  degrees of freedom
Residual deviance: 681.96  on 209  degrees of freedom
AIC: 686.82

Number of Fisher Scoring iterations: 6

valid_p_w <- predict(logit.model.w, newdata = valid, type = "response")

valid <- data.frame(valid, 'Pred_w' = valid_p_w)

head(valid_p_w)

          2           5           9          11          16          17 
0.006357329 0.391666954 0.027787765 0.019310356 0.070133111 0.010490844 

Now we can analyze our results above knowing that the predictions and model have already been adjusted for our sampling because of the weighted approach. We use the predict function to obtain the predictions from our model and make them a new vector in our validation data set.

Which approach is better? The more common approach is the weighted observation approach. Empirical simulation studies have proven that for large sample sizes (n > 1000), the weighted approach is better. In small sample sizes (n < 1000), the adjusted intercept approach is only better when you correctly specify the model. In other words, if you plan on doing variable selection because you are unsure if you have the correct variables in your model, then your model may be misspecified in its current form until after your perform variable selection. That being the case, it is probably safer to use the weighted approach. This is why most people use the weighted approach.

Stepwise

Here we will work through stepwise selection. In stepwise selection, the initial model is empty (contains no variables, only the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The most significant variable (if any) is then added to the model. Next, all variables in the model (here only one) are tested to see if they are still significant. If not, they are dropped. If so, then the remaining variables are again tested and the next most significant variable is added to the model. This process repeats until either no more significant variables are available to add or the same variable keeps being added and then removed from the model.

PROC LOGISTIC uses the selection = stepwise option to perform this technique. Notice also the slentry = and slstay = options where we specify the significance to both enter and stay in the model respectively. Here we chose 0.03 based on our sample size.

ods html select ModelBuildingSummary ParameterEstimates ORPlot;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = age race lwt smoke / selection=stepwise slentry=0.03 slstay=0.03 clodds=pl clparm=pl;
    title 'Modeling Low Birth Weight';
run;
quit;

The model building summary is shown to have only added one variable - lwt. This may be different across different stepwise techniques.

Backward

Here we will work through backward selection. In backward selection, the initial model is full (contains all variables, including the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The least significant variable (if any) is then dropped from the model. Then the remaining variables are again tested and the next least significant variable is dropped from the model. This process repeats until no more insignificant variables are available to drop.

PROC LOGISTIC uses the selection = backward option to perform this technique. Notice also the slstay = option where we specify the significance to stay in the model. Here we chose 0.03 based on our sample size.

ods html select ModelBuildingSummary ParameterEstimates ORPlot;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = age race lwt smoke / selection=backward slstay=0.03 clodds=pl clparm=pl;
    title 'Modeling Low Birth Weight';
run;
quit;

The model building summary is shown to have only dropped two variables - lwt and age. This leaves race and smoke in our model as shown in the parameter estimates table. This may be different across different stepwise techniques.

Forward with Interactions

Here we will work through forward selection. In forward selection, the initial model is empty (contains no variables, only the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The most significant variable (if any) is then added to the model. Then the remaining variables are again tested and the next most significant variable is added to the model. This process repeats until no more significant variables are available to add.

PROC LOGISTIC uses the selection = forward option to perform this technique. Notice also the slentry = option where we specify the significance to enter in the model. Here we chose 0.03 based on our sample size.

Foward selection is the least used technique because stepwise selection does the same as forward selection with the added benefit of dropping insignificant variables. The main use for forward selection is to test higher order terms and interactions in models. Here we are testing not only the variables, but their two-way interactions. This can be done in PROC LOGISTIC by placing a “|” in between variables. The “@2” is used to specify that we want only two-way interactions between these variables. The include = option forces SAS to include the first four variables - our main effects. This forces SAS to check all two-way interactions one at a time without having all of them in the model from the beginning like backward selection with interactions would do.

ods html select ModelBuildingSummary ParameterEstimates ORPlot;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = age race lwt smoke age|race|lwt|smoke @2 / selection=forward slentry=0.05 clodds=pl clparm=pl include=4;
    title 'Modeling Low Birth Weight';
run;
quit;

Notice that the modeling summary is empty as no interaction is found significant at our significance level. The model then remains at only the main effects as we told SAS to include these. From here was can use the other two techniques to select down further.

Stepwise

Here we will work through stepwise selection. In stepwise selection, the initial model is empty (contains no variables, only the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The most significant variable (if any) is then added to the model. Next, all variables in the model (here only one) are tested to see if they are still significant. If not, they are dropped. If so, then the remaining variables are again tested and the next most significant variable is added to the model. This process repeats until either no more significant variables are available to add or the same variable keeps being added and then removed from the model.

In R we must specify two models - the empty and full models - to step between. Notice the full model contains all the variables, while the empty model contains only the intercept (signified by ~ 1 in the formula). We then use the step function where we start at the empty model. The scope option tells R the lower (lower =) and upper (upper =) model to step between. The direction = "both" option tells R to use stepwise selection.

full.model <- glm(low ~ age + lwt + factor(smoke) + factor(race), 
                  data = bwt, family = binomial(link = "logit"))

empty.model <- glm(low ~ 1, data = bwt, family = binomial(link = "logit"))

step.model <- step(empty.model, 
                   scope = list(lower=formula(empty.model), 
                                upper=formula(full.model)), 
                   direction = "both")

Start:  AIC=236.67
low ~ 1

                Df Deviance    AIC
+ lwt            1   228.69 232.69
+ factor(smoke)  1   229.81 233.81
+ factor(race)   2   229.66 235.66
+ age            1   231.91 235.91
<none>               234.67 236.67

Step:  AIC=232.69
low ~ lwt

                Df Deviance    AIC
+ factor(smoke)  1   224.34 230.34
+ factor(race)   2   223.26 231.26
<none>               228.69 232.69
+ age            1   227.12 233.12
- lwt            1   234.67 236.67

Step:  AIC=230.34
low ~ lwt + factor(smoke)

                Df Deviance    AIC
+ factor(race)   2   215.01 225.01
<none>               224.34 230.34
+ age            1   222.88 230.88
- factor(smoke)  1   228.69 232.69
- lwt            1   229.81 233.81

Step:  AIC=225.01
low ~ lwt + factor(smoke) + factor(race)

                Df Deviance    AIC
<none>               215.01 225.01
+ age            1   214.58 226.58
- lwt            1   219.97 227.97
- factor(race)   2   224.34 230.34
- factor(smoke)  1   223.26 231.26

summary(step.model)

Call:
glm(formula = low ~ lwt + factor(smoke) + factor(race), family = binomial(link = "logit"), 
    data = bwt)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5278  -0.9053  -0.5863   1.2878   2.0364  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)   
(Intercept)        1.18087    1.00983   1.169  0.24225   
lwt               -0.01326    0.00631  -2.101  0.03562 * 
factor(smoke)1     1.06001    0.37832   2.802  0.00508 **
factor(race)other -0.31958    0.52560  -0.608  0.54317   
factor(race)white -1.29009    0.51087  -2.525  0.01156 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 215.01  on 184  degrees of freedom
AIC: 225.01

Number of Fisher Scoring iterations: 4

The model building summary is shown to have added three variables - lwt, smoke, and race - based on lowest AIC. This may be different across different stepwise techniques.

Backward

Here we will work through backward selection. In backward selection, the initial model is full (contains all variables, including the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The least significant variable (if any) is then dropped from the model. Then the remaining variables are again tested and the next least significant variable is dropped from the model. This process repeats until no more insignificant variables are available to drop.

In R we must specify one model - the full model - to step backward from. Notice the full model contains all the variables. We then use the step function where we start at the full model. The direction = "backward" option tells R to use backward selection.

full.model <- glm(low ~ age + lwt + factor(smoke) + factor(race), 
                  data = bwt, family = binomial(link = "logit"))

back.model <- step(full.model, direction = "backward")

Start:  AIC=226.58
low ~ age + lwt + factor(smoke) + factor(race)

                Df Deviance    AIC
- age            1   215.01 225.01
<none>               214.58 226.58
- lwt            1   218.86 228.86
- factor(race)   2   222.88 230.88
- factor(smoke)  1   222.66 232.66

Step:  AIC=225.01
low ~ lwt + factor(smoke) + factor(race)

                Df Deviance    AIC
<none>               215.01 225.01
- lwt            1   219.97 227.97
- factor(race)   2   224.34 230.34
- factor(smoke)  1   223.26 231.26

summary(back.model)

Call:
glm(formula = low ~ lwt + factor(smoke) + factor(race), family = binomial(link = "logit"), 
    data = bwt)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5278  -0.9053  -0.5863   1.2878   2.0364  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)   
(Intercept)        1.18087    1.00983   1.169  0.24225   
lwt               -0.01326    0.00631  -2.101  0.03562 * 
factor(smoke)1     1.06001    0.37832   2.802  0.00508 **
factor(race)other -0.31958    0.52560  -0.608  0.54317   
factor(race)white -1.29009    0.51087  -2.525  0.01156 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 215.01  on 184  degrees of freedom
AIC: 225.01

Number of Fisher Scoring iterations: 4

The model building summary is shown to have added three variables - lwt, smoke, and race - based on lowest AIC. This may be different across different stepwise techniques.

Forward with Interactions

Here we will work through forward selection. In forward selection, the initial model is empty (contains no variables, only the intercept). Each variable is then tried to see if it is significant based on p-value and a specified significance level. The most significant variable (if any) is then added to the model. Then the remaining variables are again tested and the next most significant variable is added to the model. This process repeats until no more significant variables are available to add.

In R we must specify two models to step between. We then use the step function where we start at the empty model. The scope option tells R the lower (lower =) and upper (upper =) model to step between. The direction = "forward" option tells R to use stepwise selection.

Foward selection is the least used technique because stepwise selection does the same as forward selection with the added benefit of dropping insignificant variables. The main use for forward selection is to test higher order terms and interactions in models. Here we are testing not only the variables, but their two-way interactions. This can be done in R by creating two models - the main effects model and the two-way interaction model. Here we set the starting point for forward selection at the main effects model and step up to the interaction model.

main.model <- glm(low ~ age + lwt + factor(smoke) + factor(race), 
                  data = bwt, family = binomial(link = "logit"))

int.model <- glm(low ~ age + lwt + factor(smoke) + factor(race) +
                       age*lwt + age*factor(smoke) + age*factor(race) +
                       lwt*factor(smoke) + lwt*factor(race) + factor(smoke)*factor(race), 
                       data = bwt, family = binomial(link = "logit"))

for.model <- step(main.model, 
                  scope = list(lower=formula(main.model), 
                               upper=formula(int.model)), 
                  direction = "forward")

Start:  AIC=226.58
low ~ age + lwt + factor(smoke) + factor(race)

                             Df Deviance    AIC
<none>                            214.58 226.58
+ lwt:factor(smoke)           1   213.66 227.66
+ age:factor(smoke)           1   214.25 228.25
+ factor(smoke):factor(race)  2   212.29 228.29
+ age:lwt                     1   214.55 228.55
+ lwt:factor(race)            2   213.18 229.18
+ age:factor(race)            2   214.00 230.00

summary(for.model)

Call:
glm(formula = low ~ age + lwt + factor(smoke) + factor(race), 
    family = binomial(link = "logit"), data = bwt)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5173  -0.9065  -0.5865   1.3035   2.0401  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)   
(Intercept)        1.564123   1.167005   1.340  0.18015   
age               -0.022478   0.034170  -0.658  0.51065   
lwt               -0.012526   0.006386  -1.961  0.04982 * 
factor(smoke)1     1.054439   0.380000   2.775  0.00552 **
factor(race)other -0.288409   0.526756  -0.548  0.58402   
factor(race)white -1.231671   0.517152  -2.382  0.01724 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 214.58  on 183  degrees of freedom
AIC: 226.58

Number of Fisher Scoring iterations: 4

Notice that the modeling summary contains no addition of any interactions as they did not reduce the AIC. The model then remains at only the main effects as we told R to include these. From here was can use the other two techniques to select down further.

SAS

PROC LOGISTIC doesn’t provide the coefficient of discrimination through an option. However, we can easily calculate it ourselves from the predicted probabilities. Using the OUTPUT statement, we gather the predicted probabilities for our training set as we have previously seen. PROC SORT is then used to sort these predicted probabilites for graphing purporses. Next, we use PROC TTEST. We aren’t actually interested in the t-test, but really the calculation of the mean difference between these two categories and the visual representation of the two histograms - the predicted probabilities of the events and non-events. In PROC TTEST we use the preddprobs data set from PROC LOGISTIC. The CLASS statement defines the two groups we are trying to compare which is our target variable. The VAR statement defines the variable we are averaging and comparing, which is our predicted probabilities phat.

proc logistic data=logistic.lowbwt noprint;
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / clodds=pl clparm=pl;
    output out=predprobs p=phat;
run;
 
proc sort data=predprobs;
    by descending low;
run;

proc ttest data=predprobs order=data;
    ods select statistics summarypanel;
    class low;
    var phat;
    title 'Coefficient of Discrimination and Plots';
run;

The coefficient of discrimination is the difference of the average which is reported on the bottom two lines of the table above. The standard error calculation is different, but our only focus is the actual mean of 0.095. We can also see the two histograms of our predictions. These histograms seems to have a lot of overlap, leading us to believe that our model has a lot events and non-events with similar predicted probabilities.

R

It is easy to calculate the difference of the average of two vectors in R which is exactly what we are doing. First, we use the predict function to gather the predictions from our model (type = "response" to get probabilities) and put them as a new column p_hat in our data frame. From there we just isolate out the events and non-events and save them in two vectors p1 and p0. This is done by conditioning the data frame to only have rows where either low == 1 or low == 0. Lastly, the coefficient of discrimination is just the difference of the mean of these two vectors. A plot fo the two histograms overlayed on each other is also provided for visual inspection.

bwt$p_hat <- predict(logit.model, type = "response")

p1 <- bwt$p_hat[bwt$low == 1]
p0 <- bwt$p_hat[bwt$low == 0]
coef_discrim <- mean(p1) - mean(p0)

ggplot(bwt, aes(p_hat, fill = factor(low))) +
  geom_density(alpha = 0.7) +
  scale_fill_grey() +
  labs(x = "Predicted Probability",
       fill = "Outcome",
       title = paste("Coefficient of Discrimination = ",
                     round(coef_discrim, 3), sep = ""))

The coefficient of discrimination is the difference of the average which is reported on the bottom two lines of the table above. The standard error calculation is different, but our only focus is the actual mean of 0.095. We can also see the two histograms of our predictions. These histograms seems to have a lot of overlap, leading us to believe that our model has a lot events and non-events with similar predicted probabilities.

SAS

SAS provides all these metrics by default in the associations table part of the PROC LOGISTIC output.

ods html select Association ParameterEstimates;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / clodds=pl clparm=pl;
    title 'Modeling Low Birth Weight';
run;
quit;

As we can see from the output, our model assigned the higher predicted probability to the observation with the low brith weight 68.3% of the time.

R

Although not provided immediately in the model summary, R can easily provide the metrics above with the Concordance and somersD functions.

Concordance(bwt$low, bwt$p_hat)

$Concordance
[1] 0.6831812

$Discordance
[1] 0.3168188

$Tied
[1] 5.551115e-17

$Pairs
[1] 7670

somersD(bwt$low, bwt$p_hat)

[1] 0.3663625

As we can see from the output, our model assigned the higher predicted probability to the observation with the low brith weight 68.3% of the time.

SAS

PROC LOGISTIC produces classification tables for any cut-off that you determine. Remember that a cut-off is simply the point where you decide to predict an event and non-event from your predicted probabilities. The ctable option in the MODEL statement creates the classification table for us. The pprob option defines the cut-off (or range of cut-offs) for this classification table. We want to look at all classification tables for cut-offs ranging from a predicted probability of 0 to 0.98 by steps of 0.02. Next, we use a DATA STEP to take the classification table from PROC LOGISTIC and calculate the Youden J statistic for each of the cut-offs from the given sensitivity and specificity provided by SAS. From there we use PROC SORT to rank these by this youden variable and print the first 5 observations with PROC PRINT.

ods html select Classification;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / ctable pprob = 0 to 0.98 by 0.02;
    ods output classification=classtable;
    title 'Modeling Low Birth Weight';
run;
quit;

data classtable;
    set classtable;
    youden = sensitivity + specificity - 100;
    drop PPV NPV Correct;
run;

proc sort data=classtable;
    by descending youden;
run;

proc print data=classtable (obs = 5);
run;

We can see that the highest Youden J statistic had a value of 24.4459. This took place at a cut-off of 0.2. Therefore, according to the Youden Index at least, the optimal cut-off for our model is 0.2. In other words, if our model predicts a probability above 0.2 then we should call this an event. Any predicted probability below 0.2 should be called a non-event.

Another commonly used visual for the balance of sensitivity and specificity across all of the cut-offs is the Receiver Operator Characteristic curve. Commonly known as the ROC curve. The ROC curve plots the balance of sensitivity vs. specificity. The “best” ROC curve is the one that reaches to the upper left hand side of the plot as that would imply that our model has both high levels of sensitivity and specificity. The worst ROC curve is represented by the diagonal line in the plot since that would imply our model is as good as randomly assigning events and non-events to our observations. This leads some to calculate the area under the ROC curve (typically called AUC) as a metric summarizing the curve itself. The math won’t be shown here, but the AUC equal to the \(c\)-statistic in the Rank-order statistics section. Isn’t math fun!?!?

SAS easily produces ROC curves from PROC LOGISTIC using the plots=ROC option.

ods html select ROCCurve;
proc logistic data=logistic.lowbwt plots(only)=ROC;
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / clodds=pl clparm=pl;
    title 'Modeling Low Birth Weight';
run;
quit;

SAS can also compare the ROC curves across many models statistically to see if they are different based on their AUC. This is done through ROC statements. The only downside is that SAS can only compare nested versions of the model in the MODEL statement. Notice how the two models in the ROC statements are just nested versions of the model in the MODEL statement. The ROCCONTRAST statement is then used with the estimate = allpairs option to statistically compare these models.

ods html select ROCContrastEstimate ROCOverlay ROCContrastTest;
proc logistic data=logistic.lowbwt plots(only)=ROC;
    class race(ref='white') / param=ref; 
    model low(event='1') = age race lwt smoke / clodds=pl clparm=pl;
    ROC 'Omit Age' race lwt smoke;
    ROC 'Omit Race' lwt smoke;
    ROCcontrast / estimate = allpairs;
    title 'Comparing ROC Curves';
run;
quit;

Although the ROC curves in the plot above look different then each other visually, statistically we cannot find any differences between them based on the p-values from the tests in the table above.

R

R produces classification tables for any cut-off that you determine. Remember that a cut-off is simply the point where you decide to predict an event and non-event from your predicted probabilities. The confusionMatrix function creates the classification table for us. The threshold option defines the cut-off for this classification table.

confusionMatrix(bwt$low, bwt$p_hat, threshold = 0.5)

    0  1
0 122 50
1   8  9

We want to look at all classification tables for cut-offs ranging from a predicted probability of 0 to 0.98 by steps of 0.02. We can easily loop through this calculation with a for loop in R. In this loop we are calculating the sensivitiy and specificity at each threshold using the sensitivity and specificity functions. We also calculate the Youden J statistic for each of the cut-offs from the given sensitivity and specificity provided by R. From there we combine these variables into a single data frame print the observations with the print function.

sens <- NULL
spec <- NULL
youden <- NULL
cutoff <- NULL

for(i in 1:49){
  cutoff = c(cutoff, i/50)
  sens <- c(sens, sensitivity(bwt$low, bwt$p_hat, threshold = i/50))
  spec <- c(spec, specificity(bwt$low, bwt$p_hat, threshold = i/50))
  youden <- c(youden, youdensIndex(bwt$low, bwt$p_hat, threshold = i/50))
}

ctable <- data.frame(cutoff, sens, spec, youden)

print(ctable)

   cutoff       sens        spec      youden
1    0.02 1.00000000 0.000000000 0.000000000
2    0.04 1.00000000 0.000000000 0.000000000
3    0.06 1.00000000 0.007692308 0.007692308
4    0.08 1.00000000 0.030769231 0.030769231
5    0.10 1.00000000 0.084615385 0.084615385
6    0.12 1.00000000 0.123076923 0.123076923
7    0.14 0.96610169 0.230769231 0.196870926
8    0.16 0.96610169 0.284615385 0.250717080
9    0.18 0.93220339 0.330769231 0.262972621
10   0.20 0.89830508 0.353846154 0.252151239
11   0.22 0.89830508 0.369230769 0.267535854
12   0.24 0.86440678 0.400000000 0.264406780
13   0.26 0.81355932 0.415384615 0.228943937
14   0.28 0.79661017 0.430769231 0.227379400
15   0.30 0.74576271 0.484615385 0.230378096
16   0.32 0.67796610 0.553846154 0.231812256
17   0.34 0.61016949 0.661538462 0.271707953
18   0.36 0.54237288 0.723076923 0.265449804
19   0.38 0.38983051 0.784615385 0.174445893
20   0.40 0.28813559 0.815384615 0.103520209
21   0.42 0.27118644 0.853846154 0.125032595
22   0.44 0.18644068 0.915384615 0.101825293
23   0.46 0.15254237 0.923076923 0.075619296
24   0.48 0.15254237 0.930769231 0.083311604
25   0.50 0.15254237 0.938461538 0.091003911
26   0.52 0.15254237 0.946153846 0.098696219
27   0.54 0.15254237 0.953846154 0.106388527
28   0.56 0.13559322 0.953846154 0.089439374
29   0.58 0.13559322 0.953846154 0.089439374
30   0.60 0.11864407 0.969230769 0.087874837
31   0.62 0.11864407 0.976923077 0.095567145
32   0.64 0.11864407 0.976923077 0.095567145
33   0.66 0.06779661 0.992307692 0.060104302
34   0.68 0.05084746 0.992307692 0.043155150
35   0.70 0.03389831 1.000000000 0.033898305
36   0.72 0.00000000 1.000000000 0.000000000
37   0.74 0.00000000 1.000000000 0.000000000
38   0.76 0.00000000 1.000000000 0.000000000
39   0.78 0.00000000 1.000000000 0.000000000
40   0.80 0.00000000 1.000000000 0.000000000
41   0.82 0.00000000 1.000000000 0.000000000
42   0.84 0.00000000 1.000000000 0.000000000
43   0.86 0.00000000 1.000000000 0.000000000
44   0.88 0.00000000 1.000000000 0.000000000
45   0.90 0.00000000 1.000000000 0.000000000
46   0.92 0.00000000 1.000000000 0.000000000
47   0.94 0.00000000 1.000000000 0.000000000
48   0.96 0.00000000 1.000000000 0.000000000
49   0.98 0.00000000 1.000000000 0.000000000

We can see that the highest Youden J statistic had a value of 0.2675. This took place at a cut-off of 0.22. Therefore, according to the Youden Index at least, the optimal cut-off for our model is 0.22. In other words, if our model predicts a probability above 0.22 then we should call this an event. Any predicted probability below 0.22 should be called a non-event.

Another commonly used visual for the balance of sensitivity and specificity across all of the cut-offs is the Receiver Operator Characteristic curve. Commonly known as the ROC curve. The ROC curve plots the balance of sensitivity vs. specificity. The “best” ROC curve is the one that reaches to the upper left hand side of the plot as that would imply that our model has both high levels of sensitivity and specificity. The worst ROC curve is represented by the diagonal line in the plot since that would imply our model is as good as randomly assigning events and non-events to our observations. This leads some to calculate the area under the ROC curve (typically called AUC) as a metric summarizing the curve itself. The math won’t be shown here, but the AUC equal to the \(c\)-statistic in the Rank-order statistics section. Isn’t math fun!?!?

R easily produces ROC curves from a variety of functions. The first combination below is the plotROC and AUROC functions that produce the ROC curve and AUC respectively.

plotROC(bwt$low, bwt$p_hat)

AUROC(bwt$low, bwt$p_hat)

[1] 0.6823338

Another common function is the performance function that produces many more plots than the ROC curve. Here the ROC curve is obtained by plotting the true positive rate by the false positive rate using the measure = "tpr" and "x.measure = "fpr" options. The AUC is also obtained from the performance function by calling the measure = "auc" option.

pred <- prediction(fitted(logit.model), factor(bwt$low))
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
plot(perf, lwd = 3, colorize = TRUE, colorkey = TRUE,
     colorize.palette = rev(gray.colors(256)))
abline(a = 0, b = 1, lty = 3)

performance(pred, measure = "auc")@y.values

[[1]]
[1] 0.685528

SAS

PROC LOGISTIC produces classification tables for any cut-off that you determine. Remember that a cut-off is simply the point where you decide to predict an event and non-event from your predicted probabilities. The ctable option in the MODEL statement creates the classification table for us. The pprob option defines the cut-off (or range of cut-offs) for this classification table. We want to look at all classification tables for cut-offs ranging from a predicted probability of 0 to 0.98 by steps of 0.02. Next, we use a DATA STEP to take the classification table from PROC LOGISTIC and calculate the F1 Score for each of the cut-offs from the given precision and recall provided by SAS. From there we use PROC SORT to rank these by this F1 variable and print the first 5 observations with PROC PRINT.

ods html select Classification;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / ctable pprob = 0 to 0.98 by 0.02;
    ods output classification=classtable;
    title 'Modeling Low Birth Weight';
run;
quit;

data classtable;
    set classtable;
    F1 = 2*(PPV*Sensitivity)/(PPV+Sensitivity);
    drop Specificity NPV Correct;
run;

proc sort data=classtable;
    by descending F1;
run;

proc print data=classtable (obs = 5);
run;

We can see that the highest F1 score had a value of 53.8071. This took place at a cut-off of 0.2. Therefore, according to the F1 score at least, the optimal cut-off for our model is 0.2. In other words, if our model predicts a probability above 0.2 then we should call this an event. Any predicted probability below 0.2 should be called a non-event.

Another common calculation using the precision of a model is the model’s lift. The lift of a model is simply calculated as the ratio of Precision to the population proportion of the event - \(Lift = PPV/\pi_1\). The interpretation of lift is really nice for explanation. Let’s imagine that your lift was 3 and your population proportion of events was 0.2. This means that in the top 20% of your customers, your model predicted 3 times the events as compared to you selecting people at random. Sometimes people plot lift charts where they plot the precision at all the different values of the population proportion (called depth).

Similar to before, we need SAS to calculate precision and recall at different cut-off values. Here we use the SCORE statement and the fitstat and outroc=roc options. This will produce the cut-offs and measures used by SAS to calculate the ROC curves. From these calculations we can derive the lift chart. Using a DATA STEP we can calculate all of the needed variables of cutoff, depth, and precision. From there we can calculate lift by using the precision as well as the population proportion of events in our low birth weight data set - here 0.3122. We then USE PROC SGPLOT to plot these.

proc logistic data=logistic.lowbwt plots(only)=(oddsratio) noprint;
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / clodds=pl clparm=pl;
    score data=logistic.lowbwt fitstat outroc=roc;
    title 'Modeling Low Birth Weight';
run;
quit;

data work.roc; 
    set work.roc; 
    cutoff = _PROB_; 
    specif = 1-_1MSPEC_; 
    depth=(_POS_+_FALPOS_)/189*100; 
    precision=_POS_/(_POS_+_FALPOS_); 
    acc=_POS_+_NEG_; 
    lift=precision/0.3122; 
run;

proc sgplot data=work.roc; 
    *where 0.005 <= depth <= 0.50; 
    series y=lift x=depth; 
    refline 1.0 / axis=y; 
    title1 "Lift Chart for Training Data"; 
    xaxis label="Depth (%)";
    yaxis label="Lift";
run; 
quit;

R

R produces classification tables for any cut-off that you determine. Remember that a cut-off is simply the point where you decide to predict an event and non-event from your predicted probabilities. The confusionMatrix function creates the classification table for us. The threshold option defines the cut-off for this classification table.

We want to look at all classification tables for cut-offs ranging from a predicted probability of 0 to 0.98 by steps of 0.02. We can easily loop through this calculation with a for loop in R. In this loop we are calculating the sensivitiy and precision at each threshold using the sensitivity and precision functions. We also calculate the F1 score for each of the cut-offs from the given sensitivity/recall and precision provided by R. From there we combine these variables into a single data frame print the observations with the print function.

prec <- NULL
reca <- NULL
f1 <- NULL
cutoff <- NULL

for(i in 1:49){
  cutoff = c(cutoff, i/50)
  reca <- c(reca, sensitivity(bwt$low, bwt$p_hat, threshold = i/50))
  prec <- c(prec, precision(bwt$low, bwt$p_hat, threshold = i/50))
  f1 <- c(f1, 2*((prec[i]*reca[i])/(prec[i]+reca[i])))
}

ctable <- data.frame(cutoff, reca, prec, f1)

print(ctable)

   cutoff       reca      prec         f1
1    0.02 1.00000000 0.3121693 0.47580645
2    0.04 1.00000000 0.3121693 0.47580645
3    0.06 1.00000000 0.3138298 0.47773279
4    0.08 1.00000000 0.3189189 0.48360656
5    0.10 1.00000000 0.3314607 0.49789030
6    0.12 1.00000000 0.3410405 0.50862069
7    0.14 0.96610169 0.3630573 0.52777778
8    0.16 0.96610169 0.3800000 0.54545455
9    0.18 0.93220339 0.3873239 0.54726368
10   0.20 0.89830508 0.3868613 0.54081633
11   0.22 0.89830508 0.3925926 0.54639175
12   0.24 0.86440678 0.3953488 0.54255319
13   0.26 0.81355932 0.3870968 0.52459016
14   0.28 0.79661017 0.3884298 0.52222222
15   0.30 0.74576271 0.3963964 0.51764706
16   0.32 0.67796610 0.4081633 0.50955414
17   0.34 0.61016949 0.4500000 0.51798561
18   0.36 0.54237288 0.4705882 0.50393701
19   0.38 0.38983051 0.4509804 0.41818182
20   0.40 0.28813559 0.4146341 0.34000000
21   0.42 0.27118644 0.4571429 0.34042553
22   0.44 0.18644068 0.5000000 0.27160494
23   0.46 0.15254237 0.4736842 0.23076923
24   0.48 0.15254237 0.5000000 0.23376623
25   0.50 0.15254237 0.5294118 0.23684211
26   0.52 0.15254237 0.5625000 0.24000000
27   0.54 0.15254237 0.6000000 0.24324324
28   0.56 0.13559322 0.5714286 0.21917808
29   0.58 0.13559322 0.5714286 0.21917808
30   0.60 0.11864407 0.6363636 0.20000000
31   0.62 0.11864407 0.7000000 0.20289855
32   0.64 0.11864407 0.7000000 0.20289855
33   0.66 0.06779661 0.8000000 0.12500000
34   0.68 0.05084746 0.7500000 0.09523810
35   0.70 0.03389831 1.0000000 0.06557377
36   0.72 0.00000000       NaN        NaN
37   0.74 0.00000000       NaN        NaN
38   0.76 0.00000000       NaN        NaN
39   0.78 0.00000000       NaN        NaN
40   0.80 0.00000000       NaN        NaN
41   0.82 0.00000000       NaN        NaN
42   0.84 0.00000000       NaN        NaN
43   0.86 0.00000000       NaN        NaN
44   0.88 0.00000000       NaN        NaN
45   0.90 0.00000000       NaN        NaN
46   0.92 0.00000000       NaN        NaN
47   0.94 0.00000000       NaN        NaN
48   0.96 0.00000000       NaN        NaN
49   0.98 0.00000000       NaN        NaN

We can see that the highest F1 score had a value of 0.5473 This took place at a cut-off of 0.18. Therefore, according to the F1 score at least, the optimal cut-off for our model is 0.18. In other words, if our model predicts a probability above 0.18 then we should call this an event. Any predicted probability below 0.18 should be called a non-event.

Another common calculation using the precision of a model is the model’s lift. The lift of a model is simply calculated as the ratio of Precision to the population proportion of the event - \(Lift = PPV/\pi_1\). The interpretation of lift is really nice for explanation. Let’s imagine that your lift was 3 and your population proportion of events was 0.2. This means that in the top 20% of your customers, your model predicted 3 times the events as compared to you selecting people at random. Sometimes people plot lift charts where they plot the precision at all the different values of the population proportion (called depth).

The performance function in R can easily calculate and plot the lift chart for us using the measure = "lift" and x.measure = "rpp" options. This plots the lift vs. the rate of positive predictions.

perf <- performance(pred, measure = "lift", x.measure = "rpp")
plot(perf, lwd = 3, colorize = TRUE, colorkey = TRUE,
     colorize.palette = rev(gray.colors(256)),
     main = "Lift Chart for Training Data")
abline(h = 1, lty = 3)

SAS

The KS statistic is not limited only two logistic regression as it can estimate the difference between any two cumulative distribution functions. Therefore, we will need an additional procedure beyond PROC LOGISTIC. First, we must obtain the predicted probabilities of our training set through the OUTPUT statement of PROC LOGISTIC. We then input this data set of predictions into PROC NPAR1WAY. With the d option we can calculate the KS D statistic. The plot = edfplot option allows us to see the cumulative distributions that are being calculated. The CLASS statement defines the variable that specifies our two groups - here the target variable low. The VAR statement defines the distribution calculation variable - here the predicted probabilities phat.

ods html select ParameterEstimates;
proc logistic data=logistic.lowbwt;
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / clodds=pl clparm=pl;
    output out=predprobs p=phat;
run;

ods html select KSTest KS2Stats EDFPlot;
proc npar1way data=predprobs d plot=edfplot;
    class low;
    var phat;
run;

From the output we can see the KS D statistic at 0.2902. The predicted probability that this occurs at (the optimal cut-off) is defined at 0.339202. In other words, if our model predicts a probability above 0.34 then we should call this an event. Any predicted probability below 0.34 should be called a non-event, according to the KS statistic.

R

The performance function in R again helps us calculate the needed metric. Using the measure = "tpr" and x.measure = "fpr" options, we can calculate the true positive rate and false positive rate across all our predictions. From there we can just use the max function to calculate the value of maximum difference between the two - the KS statistic. Finding the cut-off at this point is a little trickier with some of the needed R functions, but we essentially search the for alpha values (here the cut-offs) for the point where the KS statistic is maximized.

perf <- performance(pred, measure = "tpr", x.measure = "fpr")
KS <- max(perf@y.values[[1]] - perf@x.values[[1]])
cutoffAtKS <- unlist(perf@alpha.values)[which.max(perf@y.values[[1]] - perf@x.values[[1]])]
print(c(KS, cutoffAtKS))

[1] 0.2902216 0.3399442

plot(x = unlist(perf@alpha.values), y = (1-unlist(perf@y.values)),
     type = "l", main = "K-S Plot (EDF)",
     xlab = 'Cut-off',
     ylab = "Proportion",
     col = "red")
lines(x = unlist(perf@alpha.values), y = (1-unlist(perf@x.values)), col = "blue")

From the output we can see the KS D statistic at 0.2902. The predicted probability that this occurs at (the optimal cut-off) is defined at 0.339202. In other words, if our model predicts a probability above 0.34 then we should call this an event. Any predicted probability below 0.34 should be called a non-event, according to the KS statistic.

SAS

PROC LOGISTIC produces classification tables for any cut-off that you determine. Remember that a cut-off is simply the point where you decide to predict an event and non-event from your predicted probabilities. The ctable option in the MODEL statement creates the classification table for us. The pprob option defines the cut-off (or range of cut-offs) for this classification table. We want to look at all classification tables for cut-offs ranging from a predicted probability of 0 to 0.98 by steps of 0.02. From there we use PROC SORT to rank these by the Correct variable produced in PROC LOGISTIC and print the first 5 observations with PROC PRINT.

ods html select Classification;
proc logistic data=logistic.lowbwt plots(only)=(oddsratio);
    class race(ref='white') / param=ref; 
    model low(event='1') = race lwt smoke / ctable pprob = 0 to 0.98 by 0.02;
    ods output classification=classtable;
    title 'Modeling Low Birth Weight';
run;
quit;

proc sort data=classtable;
    by descending Correct;
run;

proc print data=classtable (obs = 5);
run;